__author__ = 'Stephanie Juneau <stephanie.juneau@noirlab.edu>, David Herrera <david.herrera@noirlab.edu>, and the Astro Data Lab Team <datalab@noirlab.edu>'
__version__ = '20230616' # yyyymmdd
__datasets__ = ['ls_dr10','sdss_dr17']
__keywords__ = ['extragalactic','galaxies','joint query','spectroscopic redshift','3d plot']
by Stéphanie Juneau, David Herrera, and the Astro Data Lab Team
In this Notebook, we explore large-scale structures of galaxies by combining spectroscopic redshifts from SDSS/BOSS with photometry from the DESI pre-imaging Legacy Survey (LS). The advantage of spectroscopic redshifts is that they are far more accurate than photometric redshifts to probe distances to galaxies (though still need to be corrected for possible distortion effects such as the finger-of-God effect, which we ignore here). The advantage of the LS photometry is that it reaches deeper than SDSS by about 1 magnitude, which yields better image quality to measure magnitudes, colors, and galaxy shapes. While there are several possible extensions to the example work included below, we will show that a simple figure of galaxy spatial locations color-coded by galaxy morphological type reveals the known morphology-density relation.
We wanted to extend indeed a little further and be able to visualize and even interact with a representation of these galaxies in the actual space. For that, we developed a 3D plot based directly (as a flat cube, without any projection or correction) on RA, DEC and z, and that is at the end of this NB.
On a technical point of view, this short notebook illustrates an example joint query between the LS DR10 photometry Tractor table, and the SDSS/BOSS DR17 specObj spectroscopy table. It uses a pre-crossmatched table based on the closest match within a 1.5 arcsec search radius.
The columns from the LS table used (Tractor) can be seen here: https://datalab.noirlab.edu/query.php?name=ls_dr10.tractor
The columns from the pre-crossmatched table can be seen here: https://datalab.noirlab.edu/query.php?name=ls_dr10.x1p5__tractor__sdss_dr17__specobj
The columns from the SDSS DR17 used for the 3D plot can be found here https://datalab.noirlab.edu/query.php?name=sdss_dr17.specobj
If you use this notebook for your published science, please acknowledge the following:
Data Lab concept paper: Fitzpatrick et al., "The NOAO Data Laboratory: a conceptual overview", SPIE, 9149, 2014, http://dx.doi.org/10.1117/12.2057445
Data Lab disclaimer: https://datalab.noirlab.edu/disclaimers.php
# std lib
from getpass import getpass
# 3rd party
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from scipy.stats import binned_statistic_2d
%matplotlib inline
from astropy.table import Table
from astropy.cosmology import Planck18 as cosmo
import plotly
import plotly.graph_objs as go
import pandas as pd
plotly.offline.init_notebook_mode()
# Data Lab
from dl import queryClient as qc
from dl import authClient as ac
Much of the functionality of Data Lab can be accessed without explicitly logging in (the service then uses an anonymous login). But some capacities, for instance saving the results of your queries to your virtual storage space, require a login (i.e. you will need a registered user account).
If you need to log in to Data Lab, un-comment the cell below and execute it:
#token = ac.login(input("Enter user name: (+ENTER) "),getpass("Enter password: (+ENTER) "))
#ac.whoAmI()
The photometry is derived from Tractor modeling of sources, and the database includes model photometry, type (shape), as well as other quantities.
The Legacy Survey DR10 database is called ls_dr10 and includes several tables. We will use the tractor table together with the positional crossmatched table with specObj table from SDSS/BOSS DR17. The column names and descriptions can be found from the Data Lab Query Interface or using the Table Access Protocol (TAP) service with, e.g., TOPCAT.
The SDSS DR17 database is called sdss_dr17 and also includes several tables. We will use the specobj table, which has spectroscopic information.
%%time
# number of rows from LS DR3 tractor (NOTE: tractor is the main photometry table):
query="SELECT nrows FROM tbl_stat WHERE schema='ls_dr10' and tbl_name='tractor'"
# Call query manager
response = qc.query(sql=query, fmt='csv')
print(response)
nrows 3145841920 CPU times: user 25.3 ms, sys: 2.82 ms, total: 28.1 ms Wall time: 145 ms
%%time
# number of rows from SDSS specObj DR13:
query="SELECT nrows FROM tbl_stat WHERE schema='sdss_dr16' and tbl_name='specobj'"
# Call query manager
response = qc.query(sql=query, fmt='csv')
print(response)
nrows 5107041 CPU times: user 21.2 ms, sys: 1.75 ms, total: 23 ms Wall time: 83.8 ms
# ls_dr10.tractor # LS DR10 tractor photometry
# ls_dr10.x1p5__tractor__sdss_dr17__specobj # LS DR10 pre-crossmatched to SDSS DR17 specobj
# sdss_dr17.specobj # SDSS DR17 specobj
# Write query statement (sql)
query = ("""
SELECT L.ra, L.dec, L.type, L.sersic, L.g_r, L.r_z,
S.z, S.ra as plug_ra, S.dec as plug_dec, S.class, S.veldisp, S.veldisperr
FROM ls_dr10.tractor AS L
JOIN ls_dr10.x1p5__tractor__sdss_dr17__specobj AS X ON L.ls_id=X.id1
JOIN sdss_dr17.specobj AS S ON X.id2 = S.specobjid
WHERE L.ra BETWEEN %s AND %s AND L.dec BETWEEN %s AND %s AND (L.ra_ivar > 0)
LIMIT 100000
""") % (126,131,7.,12.) # small region
# L.ra, L.dec = RA, Dec from Legacy Survey (LS) table
# L.type = object type (PSF, SIMP, EXP, DEV, COMP)
# L.g_r, L.r_z = pre-computed g-r and r-z colors from photometry
# S.z = redshift (z) from SDSS specObj table
# S.plug_ra, dec = RA, Dec of SDSS fiber from specObj table
# S.class = Source class (Star, Galaxy, QSO) from SDSS
# S.veldisp, veldisp = velocity dispersion (and error) from SDSS specObj table
#
# WHERE: requirement that RA & Dec coordinates are within a rectangular region
print(query)
SELECT L.ra, L.dec, L.type, L.sersic, L.g_r, L.r_z,
S.z, S.ra as plug_ra, S.dec as plug_dec, S.class, S.veldisp, S.veldisperr
FROM ls_dr10.tractor AS L
JOIN ls_dr10.x1p5__tractor__sdss_dr17__specobj AS X ON L.ls_id=X.id1
JOIN sdss_dr17.specobj AS S ON X.id2 = S.specobjid
WHERE L.ra BETWEEN 126 AND 131 AND L.dec BETWEEN 7.0 AND 12.0 AND (L.ra_ivar > 0)
LIMIT 100000
%%time
# Call query client and save output as Astropy Table
result = qc.query(adql=query, fmt='table')
CPU times: user 86.7 ms, sys: 13.7 ms, total: 100 ms Wall time: 20.7 s
# Print length of table (Nb of rows) and the first 10 rows
print(len(result))
result[:10]
8555
| ra | dec | type | sersic | g_r | r_z | z | plug_ra | plug_dec | class | veldisp | veldisperr |
|---|---|---|---|---|---|---|---|---|---|---|---|
| float64 | float64 | str3 | float64 | float64 | float64 | float64 | float64 | float64 | str6 | float64 | float64 |
| 126.0000864493582 | 7.034927726600151 | SER | 6.0 | 1.3086281 | 0.7673378 | 0.21285242 | 126.00008000000003 | 7.0349361 | GALAXY | 276.82166 | 15.253322 |
| 126.0149353747927 | 7.007741109416743 | SER | 5.506159 | 1.3768082 | 0.77913666 | 0.24393715 | 126.01495 | 7.0077364 | GALAXY | 195.84337 | 17.137558 |
| 126.0266308179383 | 7.045709718170351 | SER | 6.0 | 1.4943676 | 0.78222275 | 0.2591496 | 126.02663000000001 | 7.0457033 | GALAXY | 218.10497 | 12.376819 |
| 126.0241703658048 | 7.067103011297897 | SER | 0.9314089 | 1.1323223 | 0.9355583 | 0.06391467 | 126.02422 | 7.0671055 | GALAXY | 144.77763 | 10.631522 |
| 126.01353786704 | 7.07235585789409 | SER | 2.2828588 | 0.9481163 | 0.70143414 | 0.064374804 | 126.01356 | 7.0723868 | GALAXY | 212.39644 | 10.097903 |
| 126.0035506525818 | 7.103901614620744 | PSF | 0.0 | 0.89271545 | 0.6074829 | 1.2897644 | 126.00357 | 7.1038999 | QSO | 0.0 | 0.0 |
| 126.0837762847585 | 7.062631211176722 | SER | 6.0 | 1.5833492 | 0.8626728 | 0.31347042 | 126.08379 | 7.0626316 | GALAXY | 273.91946 | 22.274921 |
| 126.1102945352017 | 7.059559957377497 | PSF | 0.0 | 0.013511658 | 0.344162 | 2.2276993 | 126.11029000000002 | 7.0595677 | QSO | 0.0 | 0.0 |
| 126.1333879010724 | 7.042289793306823 | PSF | 0.0 | 0.28848648 | 0.035743713 | 0.0002704559 | 126.13341 | 7.0422755 | STAR | 0.0 | 0.0 |
| 126.1386038058223 | 7.043977744211515 | SER | 2.5463092 | 1.94207 | 1.1849613 | 0.50945634 | 126.13862 | 7.0439339 | GALAXY | 339.53528 | 66.76817 |
# convert RA coordinates from [0,360] to [-180,180]
chgsign = np.where(result['ra'] > 180)
result['ra'][chgsign] = result['ra'][chgsign]-360.
result['plug_ra'][chgsign] = result['plug_ra'][chgsign]-360.
plt.figure(figsize=(9,8))
# plot RA, Dec from LS catalog in red with larger symbols
plt.scatter(result['ra'],result['dec'],s=3,color='red',marker='1')
# overplot RA, Dec from SDSS catalog in blue with smaller symbols
plt.scatter(result['plug_ra'],result['plug_dec'],s=3,color='black',marker='2')
# Extent of RA, Dec (in degrees) to plot
xmin = 126.
xmax = 131.
ymin = 7.
ymax = 12.
plt.axis([xmin, xmax, ymin, ymax])
plt.xlim(reversed(plt.xlim())) # flip the x-axis
plt.xlabel("RA (degrees)", fontsize=20)
plt.ylabel("Dec (degrees)", fontsize=20)
plt.show()
Plot the positions of a broad range of redshift, and overplot a thin slice in redshift to show possible structures within that slice.
# Select redshift slice
rz = (result['z'] >0.105) & (result['z']<0.125)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6.5))
# plot all points in red (all redshifts)
ax1.scatter(result['plug_ra'],result['plug_dec'],s=3,color='r',marker='o',alpha=0.25)
# overplot in blue objects in narrow redshift slice
ax1.scatter(result['plug_ra'][rz],result['plug_dec'][rz],s=10,color='black')
# Extent of RA, Dec (in degrees) to plot
xmin = 126.
xmax = 131.
ymin = 7.
ymax = 12.
ax1.axis([xmin, xmax, ymin, ymax])
ax1.set_xlim(reversed(ax1.set_xlim())) # flip the x-axis
ax1.set_xlabel("RA (degrees)", fontsize=20)
ax1.set_ylabel("Dec (degrees)", fontsize=20)
# add rectangle to show where we zoom in next panel
ax1.add_patch(patches.Rectangle((128.65-0.25, 8.85-0.2),0.5,0.4,fill=False,color='b'))
## ZOOM IN A SMALLER REGION
# plot all points in red (all redshifts)
ax2.scatter(result['plug_ra'],result['plug_dec'],s=15,color='r',marker='o',alpha=0.3)
# overplot in blue objects in narrow redshift slice
ax2.scatter(result['plug_ra'][rz],result['plug_dec'][rz],s=30,color='black')
# Extent of RA, Dec (in degrees) to plot
xmin = 128.4
xmax = 128.9
ymin = 8.65
ymax = 9.05
ax2.axis([xmin, xmax, ymin, ymax])
ax2.set_xlim(reversed(ax2.set_xlim())) # flip the x-axis
ax2.set_xlabel("RA (degrees)", fontsize=20)
ax2.set_ylabel("Dec (degrees)", fontsize=20)
# add rectangle to show where we zoom in next panel
ax2.add_patch(patches.Rectangle((128.63, 8.81),0.13,0.107,fill=False,color='b'))
plt.show()
Above, the left-hand panel shows a thin redshift slice (0.105 < z < 0.125, black symbols) among objects with redshifts from the spectroscopic SDSS DR13 sample (red symbols). We can see by eye some large-scale filamentary structures and overdensities. The blue rectangle shows a selected region where we zoom in the right-hand panel. On the latter, we further select a smaller region, which we will use in the next cell below.
## ZOOM IN AGAIN OVER AN EVEN SMALLER REGION
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5.5))
# plot all points in red (all redshifts)
ax1.scatter(result['plug_ra'],result['plug_dec'],s=30,color='r',marker='o',alpha=0.25)
# overplot in black objects in narrow redshift slice
ax1.scatter(result['plug_ra'][rz],result['plug_dec'][rz],s=50,color='black')
# Extent of RA, Dec (in degrees) to plot
xmin = 128.63
xmax = 128.76
ymin = 8.81
ymax = 8.917
ax1.axis([xmin, xmax, ymin, ymax])
ax1.set_xlim(reversed(ax1.set_xlim())) # flip the x-axis
ax1.set_xlabel("RA (degrees)", fontsize=20)
ax1.set_ylabel("Dec (degrees)", fontsize=20)
## SHOW DECaLS IMAGE (screenshot pre-made but could instead implement image cutout)
im = plt.imread('DECaLS_screenshot_zoomIn_labels.jpg')
ax2.imshow(im)
ax2.axis('off')
plt.show()
The left-hand panel shows the small region enclosed in the blue rectangle that we chose above (right-hand panel). The galaxies in black are in the same narrow redshift slice as defined previously (0.105 < z < 0.125). The right-hand panel is an image cutout of the same region of the sky from the LS sky viewer. The galaxies encircled correspond to the points in black, and some or perhaps most of them likely belong to a galaxy cluster.
There are many possible extensions to this work. For instance, one could plot again with symbols coded with object type (from LS) and/or class (from SDSS) and/or velocity dispersion (from SDSS) and/or other quantities. Here, we will start with the object "TYPE" from LS, related to the morphological shapes.
The object shape (2D light profile) is modeled by the Tractor (Lang, Hogg & Mykytyn) as part of the procedure to compute model photometry.
Possible shapes for LS DR10 Morphological Classification:
# Select redshift slice
rz = (result['z'] >0.105) & (result['z']<0.125)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6.5))
# plot all points in red (all redshifts)
ax1.scatter(result['plug_ra'],result['plug_dec'],s=3,color='r',marker='o',alpha=0.25)
# overplot in black objects in narrow redshift slice
ax1.scatter(result['plug_ra'][rz],result['plug_dec'][rz],s=10,color='black')
# Extent of RA, Dec (in degrees) to plot
xmin = 126.
xmax = 131.
ymin = 7.
ymax = 12.
ax1.axis([xmin, xmax, ymin, ymax])
ax1.set_xlim(reversed(ax1.set_xlim())) # flip the x-axis
ax1.set_xlabel("RA (degrees)", fontsize=20)
ax1.set_ylabel("Dec (degrees)", fontsize=20)
# Select redshift slice and per morphological type
# DEV or early type Sersic with steep index Sersic>=3
rdev = np.where(((result['type']=='DEV')|((result['type']=='SER')&(result['sersic']>=3))) & rz)
# EXP or late type Sersic with index Sersic<2
rexp = np.where(((result['type']=='EXP')|((result['type']=='SER')&(result['sersic']<2))) & rz)
# Intermediate with 2<= Sersic < 3
rint = np.where(((result['type']=='SER')&(result['sersic']>=2)&(result['sersic']<3)) & rz)
# plot all points in red (all redshifts)
ax2.scatter(result['plug_ra'],result['plug_dec'],s=3,color='gray',marker='o',alpha=0.25)
# overplot in objects in narrow redshift slice color-coded blue (EXP or Late type) or red (DEV or Early type)
ax2.scatter(result['plug_ra'][rexp],result['plug_dec'][rexp],s=15,color='b') # blue = EXP or late-type Sersic
ax2.scatter(result['plug_ra'][rint],result['plug_dec'][rint],s=15,color='orange') # orange = Intermediate
ax2.scatter(result['plug_ra'][rdev],result['plug_dec'][rdev],s=15,color='r') # red = DEV or early-type Sersic
plt.axis([xmin, xmax, ymin, ymax])
plt.xlim(reversed(plt.xlim())) # flip the x-axis
plt.xlabel("RA (degrees)", fontsize=20)
plt.ylabel("Dec (degrees)", fontsize=20)
plt.show()
There are pre-computed colors available. The columns are described here: https://datalab.noirlab.edu/query.php?name=ls_dr10.tractor
Another possibility would be to plot again the galaxies spatial coordinates, but color-coded according to their photometric colors. This is left as an exercise for the user, but feel free to get in touch with the Astro Data Lab Team if you have questions.
We can explore filaments and clusters of galaxies better if we can plot them in 3D. We turn to a different area of the sky, and will query for a sample of SDSS galaxies in the near to slightly distant universe, and plot in 3D the cone containing them. To avoid conamination by the galactic plane of the Milky Way, we point our search cone at high galactic latitudes.
We select (mostly) SDSS galaxies within a 10-degree radius around a high galactic latitude direction, (ra,dec) = (160,45) degrees. We limit our search to positive redshifts between 0.02 and 0.3.
# Create the query to fetch the SDSS data from DataLab:
query = """
SELECT ra,dec,z
FROM sdss_dr17.specobj
WHERE q3c_radial_query(ra,dec,%s,%s,%s)
AND z between 0.02 AND 0.3
""" % (160,45,10)
print (query)
SELECT ra,dec,z FROM sdss_dr17.specobj WHERE q3c_radial_query(ra,dec,160,45,10) AND z between 0.02 AND 0.3
Run the query to fetch the SDSS data from the ls_dr17.specobj table
%%time
selection = qc.query(adql=query, fmt='csv')
CPU times: user 25.6 ms, sys: 14.3 ms, total: 39.9 ms Wall time: 1.28 s
Reformat output into a table
data = Table.read(selection, format='csv') #Astropy Table
print("Number of galaxies in the sample: %d" % len(data))
data[:5]
Number of galaxies in the sample: 36331
| ra | dec | z |
|---|---|---|
| float64 | float64 | float64 |
| 159.25732 | 54.941463 | 0.072250396 |
| 159.17507999999998 | 54.942007 | 0.14156331 |
| 162.24439 | 54.889704 | 0.14798966 |
| 162.50349 | 54.858997 | 0.17074445 |
| 162.42234 | 54.829701 | 0.18014126 |
Compute luminosity distance for every galaxy, using Planck2018 cosmology values
dist = cosmo.luminosity_distance(data['z'])
dist
To plot the positions of each galaxy in true 3D space, we will convert the angular coordinates ra & dec, and the distance coordinate, to Cartesian coordinates X,Y,Z (all measured in Mpc from the coordinate system origin).
def get_cartesian(ra,dec,dist):
# convert ra and dec to radians, since numpy expects this as arguments to trigonometric functions
rarad = np.radians(ra)
decrad = np.radians(dec)
X = dist * np.sin(decrad) * np.cos(rarad)
Y = dist * np.sin(decrad) * np.sin(rarad)
Z = dist * np.cos(decrad)
return X,Y,Z
X,Y,Z = get_cartesian(data['ra'],data['dec'],dist)
To create an interactive 3D plot (one where we can zoom, pan, and rotate the scene), we will use the plotly package. The galaxies in our sample will be plotted using a 3D scatter routine. We also set some overall properties of the plot, such as the size of the markers, and a color map (we color each galaxy redder the further away it is).
trace = go.Scatter3d(
x = X,
y = Y,
z = Z,
mode = 'markers',
marker = {
'size' : 0.7,
'opacity' : 0.5,
'color' : dist,
'colorscale': 'OrRd'
}
)
data = [trace]
Next, we define camera location and the layout of the plot
# set up the view point
camera = dict(
up = dict(x = 0, y = 0, z = 1),
center = dict(x = 0.2, y = 0, z = 0),
eye = dict(x = 0.6, y = -0.6, z = -1.0)
)
# set up the plot scene
layout = go.Layout(
scene = dict(
xaxis = dict(title = 'X',
backgroundcolor = 'black',
gridcolor = "rgb(40,40,40)"),
yaxis = dict(title = 'Y',
backgroundcolor = 'black',
gridcolor = "rgb(40,40,40)"),
zaxis = dict(title = 'Z',
backgroundcolor = 'black',
gridcolor = "rgb(40,40,40)"),
),
scene_camera = camera,
plot_bgcolor = 'black',
paper_bgcolor = 'black',
title = None,
showlegend = False,
width = 800,
height = 800,
autosize = False,
margin = {'l':0, 'r':0, 'b':0, 't':0},
dragmode = 'orbit'
)
# Draw the plot
plot_figure = go.Figure(data=data, layout=layout)
plot_figure.update_layout()
Explore the scene by:
Observe how the galaxies form compact groups and clusters, and on larger scales form huge filaments. This is the large scale structure of the universe. In details, we can also notice some elongation along the direction of redshift, which is sometimes called the "Fingers of God" effect. This is due to additional contributions to the redshift coming from the proper motions of galaxies. This effect is usually corrected in a statistical manner rather than on a galaxy-per-galaxy basis (other useful explainer).